Fit temperature models and predict growing season temperature

Author

Max Lindmark, Jan Ohlberger, Anna Gårdmark

Published

September 22, 2023

Load libraries

pkgs <- c("here","tidyverse", "tidylog", "RColorBrewer", "viridis", "sdmTMB", "sdmTMBextra", "patchwork", "RCurl", "tidylog") 

# minpack.lm needed if using nlsLM()
if(length(setdiff(pkgs,rownames(installed.packages()))) > 0){

    install.packages(setdiff(pkgs, rownames(installed.packages())), dependencies = T)
  
  }

invisible(lapply(pkgs, library,character.only = T))
here() starts at /Users/maxlindmark/Dropbox/Max work/R/perch-growth
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.2     ✔ readr     2.1.4
✔ forcats   1.0.0     ✔ stringr   1.5.0
✔ ggplot2   3.4.2     ✔ tibble    3.2.1
✔ lubridate 1.9.2     ✔ tidyr     1.3.0
✔ purrr     1.0.1     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

Attaching package: 'tidylog'


The following objects are masked from 'package:dplyr':

    add_count, add_tally, anti_join, count, distinct, distinct_all,
    distinct_at, distinct_if, filter, filter_all, filter_at, filter_if,
    full_join, group_by, group_by_all, group_by_at, group_by_if,
    inner_join, left_join, mutate, mutate_all, mutate_at, mutate_if,
    relocate, rename, rename_all, rename_at, rename_if, rename_with,
    right_join, sample_frac, sample_n, select, select_all, select_at,
    select_if, semi_join, slice, slice_head, slice_max, slice_min,
    slice_sample, slice_tail, summarise, summarise_all, summarise_at,
    summarise_if, summarize, summarize_all, summarize_at, summarize_if,
    tally, top_frac, top_n, transmute, transmute_all, transmute_at,
    transmute_if, ungroup


The following objects are masked from 'package:tidyr':

    drop_na, fill, gather, pivot_longer, pivot_wider, replace_na,
    spread, uncount


The following object is masked from 'package:stats':

    filter


Loading required package: viridisLite


Attaching package: 'sdmTMBextra'


The following objects are masked from 'package:sdmTMB':

    dharma_residuals, extract_mcmc



Attaching package: 'RCurl'


The following object is masked from 'package:tidyr':

    complete
# devtools::install_github("seananderson/ggsidekick") # not on CRAN 
library(ggsidekick)
theme_set(theme_sleek())

# Set path:
home <- here::here()

Load cache

# qwraps2::lazyload_cache_dir(path = paste0(home, "/R/analyze-data/01-fit-temp-models-predict_cache/html"))

Read data

d <- read_csv(paste0(home, "/output/temp_data_for_fitting.csv"))
Rows: 98616 Columns: 13
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (4): area, source, db, id
dbl  (8): year, temp, yday, month, day, VkDag, stn_nr, section_nr
date (1): date

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
d <- d |> mutate(area = as.factor(area),
                 source_f = as.factor(source),
                 year_f = as.factor(year)) |> 
  filter(!area %in% c("VN", "TH")) # VN: no logger data, TH: to short time series
mutate: converted 'area' from character to factor (0 new NA)
        new variable 'source_f' (factor) with 3 unique values and 0% NA
        new variable 'year_f' (factor) with 83 unique values and 0% NA
filter: removed 5,545 rows (6%), 93,071 rows remaining
# Keep track of the years for which we have cohorts for matching with cohort data
gdat <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/perch-growth/master/data/for-analysis/dat.csv")
Rows: 338460 Columns: 12
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (6): age_ring, area, gear, ID, sex, keep
dbl (6): length_mm, age_bc, age_catch, catch_year, cohort, final_length

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
gdat$area_cohort_age <- as.factor(paste(gdat$area, gdat$cohort, gdat$age_bc))

gdat <- gdat |>
  group_by(area_cohort_age) |> 
  filter(n() > 10) |> 
  filter(age_catch > 3) |> 
  group_by(area) |>
  summarise(min = min(cohort),
            max = max(cohort)) |> 
  arrange(min)
group_by: one grouping variable (area_cohort_age)
filter (grouped): removed 5,190 rows (2%), 333,270 rows remaining
filter (grouped): removed 107,058 rows (32%), 226,212 rows remaining
group_by: one grouping variable (area)
summarise: now 12 rows and 3 columns, ungrouped
d <- left_join(d, gdat, by = "area") |>
  mutate(area = as.factor(area),
         growth_dat = ifelse(year >= min & year <= max, "Y", "N"))
left_join: added 2 columns (min, max)
           > rows only in x        0
           > rows only in y  (     2)
           > matched rows     93,071
           >                 ========
           > rows total       93,071
mutate: converted 'area' from character to factor (0 new NA)
        new variable 'growth_dat' (character) with 2 unique values and 0% NA
# Drop data in SI_HA and BT before onset of warming
d <- d |>
  mutate(discard = "N",
         discard = ifelse(area == "BT" & year <= 1980, "Y", discard),
         discard = ifelse(area == "SI_HA" & year <= 1972, "Y", discard)) |> 
  filter(discard == "N")
mutate: new variable 'discard' (character) with 2 unique values and 0% NA
filter: removed 1,146 rows (1%), 91,925 rows remaining
# Drop heated areas actually for the full plot
df <- d |> filter(!area %in% c("BT", "SI_HA"))
filter: removed 24,149 rows (26%), 67,776 rows remaining

Fit models

Source as independent or interactive effect

m <- sdmTMB(temp ~ area*year_f + source_f + s(yday, by = area, bs = "cc"), 
            data = df,
            family = student(df = 5),
            spatial = "off",
            spatiotemporal = "off",
            knots = list(yday = c(0.5, 364.5)),
            control = sdmTMBcontrol(newton_loops = 1))

Try alternative degrees of freedom (not evaluated)

#df=20 (effectively gaussian)
m2 <- sdmTMB(temp ~ area*year_f + source_f + s(yday, by = area, bs = "cc"), 
            data = d,
            family = student(df = 20),
            spatial = "off",
            spatiotemporal = "off",
            knots = list(yday = c(0.5, 364.5)),
            control = sdmTMBcontrol(newton_loops = 1))

# df = 9
m3 <- sdmTMB(temp ~ area*year_f + source_f + s(yday, by = area, bs = "cc"), 
            data = d,
            family = student(df = 9),
            spatial = "off",
            spatiotemporal = "off",
            knots = list(yday = c(0.5, 364.5)),
            control = sdmTMBcontrol(newton_loops = 1))

# df=4
m4 <- sdmTMB(temp ~ area*year_f + source_f + s(yday, by = area, bs = "cc"), 
            data = d,
            family = student(df = 4),
            spatial = "off",
            spatiotemporal = "off",
            knots = list(yday = c(0.5, 364.5)),
            control = sdmTMBcontrol(newton_loops = 1))

# Plot all residuals
mcmc_res <- residuals(m, type = "mle-mcmc",
                      mcmc_samples = sdmTMBextra::predict_mle_mcmc(m,
                                                                   mcmc_iter = 201,
                                                                   mcmc_warmup = 200))

mcmc_res20 <- residuals(m2, type = "mle-mcmc",
                        mcmc_samples = sdmTMBextra::predict_mle_mcmc(m2,
                                                                     mcmc_iter = 201,
                                                                     mcmc_warmup = 200))

mcmc_res9 <- residuals(m3, type = "mle-mcmc",
                       mcmc_samples = sdmTMBextra::predict_mle_mcmc(m3,
                                                                    mcmc_iter = 201,
                                                                    mcmc_warmup = 200))

mcmc_res4 <- residuals(m4, type = "mle-mcmc",
                       mcmc_samples = sdmTMBextra::predict_mle_mcmc(m4,
                                                                    mcmc_iter = 201,
                                                                    mcmc_warmup = 200))

dres <- df |> mutate("df=5" = mcmc_res,
                     "df=20" = mcmc_res20,
                     "df=9" = mcmc_res9,
                     "df=4" = mcmc_res4) |> 
  select(`df=5`, `df=20`, `df=9`, `df=4`) |> 
  pivot_longer(everything())

ggplot(dres, aes(sample = value)) +
  stat_qq(size = 0.75, shape = 21, fill = NA) +
  facet_wrap(~factor(name, levels = c("df=20", "df=9", "df=5", "df=4"))) +
  stat_qq_line() +
  labs(y = "Sample Quantiles", x = "Theoretical Quantiles") + 
  theme(aspect.ratio = 1)

Check fit

sanity(m)
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
mcmc_res <- residuals(m, type = "mle-mcmc",
                      mcmc_samples = sdmTMBextra::predict_mle_mcmc(m,
                                                                   mcmc_iter = 201,
                                                                   mcmc_warmup = 200))

SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.007709 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 77.09 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 86.1689 seconds (Warm-up)
Chain 1:                0.101675 seconds (Sampling)
Chain 1:                86.2705 seconds (Total)
Chain 1: 
ggplot(df, aes(sample = mcmc_res)) +
  stat_qq() +
  stat_qq_line() +
  labs(y = "Sample Quantiles", x = "Theoretical Quantiles") + 
  theme(aspect.ratio = 1)

ggsave(paste0(home, "/figures/supp/full_model/qq_temp.pdf"), width = 11, height = 11, units = "cm")

summary(m)
Model fit by ML ['sdmTMB']
Formula: temp ~ area * year_f + source_f + s(yday, by = area, bs = "cc")
Mesh: NULL (isotropic covariance)
Data: df
Family: student(link = 'identity')
 
                     coef.est coef.se
(Intercept)              8.83    0.47
areaFB                   0.07    0.74
areaFM                   0.13    0.78
areaHO                  -0.21    0.81
areaJM                   0.68    0.76
areaMU                   1.31    0.73
areaRA                  -0.79    0.80
areaSI_EK                1.96    0.74
year_f1941              -0.28    0.69
year_f1942              -0.29    0.67
year_f1943               1.34    0.68
year_f1944               1.35    0.68
year_f1945               1.69    0.67
year_f1946               1.92    0.70
year_f1947               0.86    0.69
year_f1948               1.40    0.67
year_f1949               1.38    0.67
year_f1950               1.16    0.66
year_f1951               1.05    0.68
year_f1952               0.22    0.67
year_f1953               1.62    0.67
year_f1954               1.62    0.67
year_f1955               1.57    0.68
year_f1956               1.07    0.68
year_f1957               1.25    0.68
year_f1958               1.14    0.69
year_f1959               2.26    0.67
year_f1960               1.27    0.66
year_f1961               2.36    0.69
year_f1962               1.42    0.70
year_f1963               1.58    0.67
year_f1964               1.33    0.68
year_f1965               1.25    0.68
year_f1966               1.08    0.66
year_f1967               1.93    0.67
year_f1968               1.63    0.66
year_f1969               1.27    0.65
year_f1970               1.25    0.66
year_f1971               1.19    0.66
year_f1972               1.43    0.67
year_f1973               2.02    0.71
year_f1974               1.87    0.66
year_f1975               2.23    0.66
year_f1976               0.60    0.66
year_f1977               0.62    0.66
year_f1978               0.47    0.66
year_f1979               0.10    0.67
year_f1980               0.57    0.67
year_f1981               0.54    0.66
year_f1982               0.73    0.67
year_f1983               1.30    0.66
year_f1984               1.29    0.66
year_f1985              -0.11    0.66
year_f1986               0.39    0.67
year_f1987              -0.82    0.67
year_f1988               1.34    0.67
year_f1989               2.10    0.67
year_f1990               2.33    0.69
year_f1991               1.55    0.54
year_f1992               1.08    0.51
year_f1993               0.84    0.51
year_f1994               1.31    0.50
year_f1995               2.26    0.52
year_f1996               0.69    0.51
year_f1997               4.63    0.52
year_f1998               0.68    0.48
year_f1999               2.28    0.49
year_f2000               2.03    0.48
year_f2001               2.79    0.48
year_f2002               4.01    0.49
year_f2003               3.28    0.50
year_f2004               2.15    0.48
year_f2005               2.82    0.49
year_f2006               0.93    0.53
year_f2007               1.35    0.52
year_f2008               2.20    0.53
year_f2009               2.52    0.53
year_f2010               1.34    0.68
year_f2011               1.75    0.69
year_f2012               1.74    0.66
year_f2013               2.13    0.67
year_f2014               2.50    0.66
year_f2015               2.55    0.66
year_f2016               2.65    0.66
year_f2017               2.20    0.66
year_f2018               2.80    0.67
year_f2019               2.14    0.66
year_f2020               3.22    0.65
year_f2021               2.34    0.67
year_f2022               2.27    1.21
source_ffishing          0.33    0.05
source_flogger           0.19    0.04
areaFB:year_f1941        0.05    1.13
areaFM:year_f1941        0.06    1.19
areaHO:year_f1941        0.41    1.20
areaJM:year_f1941       -0.19    1.16
areaMU:year_f1941       -0.16    1.09
areaRA:year_f1941        0.10    1.16
areaSI_EK:year_f1941    -0.44    1.11
areaFB:year_f1942        0.14    1.06
areaFM:year_f1942        0.15    1.11
areaHO:year_f1942        0.38    1.14
areaJM:year_f1942       -0.05    1.10
areaMU:year_f1942       -0.19    1.02
areaRA:year_f1942        0.05    1.10
areaSI_EK:year_f1942    -0.19    1.08
areaFB:year_f1943       -0.57    1.06
areaFM:year_f1943       -0.78    1.11
areaHO:year_f1943       -0.86    1.15
areaJM:year_f1943       -0.64    1.09
areaMU:year_f1943       -0.40    1.02
areaRA:year_f1943       -0.85    1.15
areaSI_EK:year_f1943    -0.31    1.04
areaFB:year_f1944       -0.32    1.09
areaFM:year_f1944       -0.31    1.16
areaHO:year_f1944       -0.43    1.19
areaJM:year_f1944       -0.41    1.12
areaMU:year_f1944       -0.24    1.03
areaRA:year_f1944       -0.78    1.15
areaSI_EK:year_f1944    -0.15    1.06
areaFB:year_f1945       -0.23    1.04
areaFM:year_f1945       -0.42    1.09
areaHO:year_f1945       -0.74    1.13
areaJM:year_f1945       -0.18    1.07
areaMU:year_f1945       -0.16    1.00
areaRA:year_f1945       -1.15    1.13
areaSI_EK:year_f1945     0.07    1.03
areaFB:year_f1946        0.13    1.05
areaFM:year_f1946       -0.07    1.11
areaHO:year_f1946       -0.58    1.14
areaJM:year_f1946        0.19    1.09
areaMU:year_f1946       -0.08    1.01
areaRA:year_f1946       -1.67    1.13
areaSI_EK:year_f1946    -0.01    1.04
areaFB:year_f1947       -0.10    1.05
areaFM:year_f1947       -0.23    1.11
areaHO:year_f1947       -1.00    1.14
areaJM:year_f1947        0.03    1.10
areaMU:year_f1947       -0.02    1.02
areaRA:year_f1947       -1.85    1.19
areaSI_EK:year_f1947     0.21    1.07
areaFB:year_f1948       -0.11    1.04
areaFM:year_f1948       -0.15    1.11
areaHO:year_f1948       -0.23    1.13
areaJM:year_f1948       -0.13    1.07
areaMU:year_f1948       -0.27    1.01
areaRA:year_f1948       -0.84    1.10
areaSI_EK:year_f1948    -0.26    1.05
areaFB:year_f1949       -0.20    1.09
areaFM:year_f1949       -0.39    1.15
areaHO:year_f1949       -0.54    1.21
areaJM:year_f1949       -0.17    1.12
areaMU:year_f1949        0.06    1.04
areaRA:year_f1949       -1.36    1.18
areaSI_EK:year_f1949     0.31    1.09
areaFB:year_f1950        0.01    1.05
areaFM:year_f1950       -0.12    1.12
areaHO:year_f1950       -0.22    1.15
areaJM:year_f1950        0.22    1.08
areaMU:year_f1950        0.01    1.00
areaRA:year_f1950       -0.60    1.15
areaSI_EK:year_f1950     0.24    1.02
areaFB:year_f1951       -0.27    1.12
areaFM:year_f1951       -0.40    1.20
areaHO:year_f1951       -0.43    1.22
areaJM:year_f1951       -0.30    1.16
areaMU:year_f1951       -0.03    1.06
areaRA:year_f1951       -1.11    1.18
areaSI_EK:year_f1951     0.12    1.11
areaFB:year_f1952       -0.11    1.03
areaFM:year_f1952       -0.23    1.07
areaHO:year_f1952       -0.40    1.15
areaJM:year_f1952        0.15    1.06
areaMU:year_f1952       -0.10    1.01
areaRA:year_f1952       -0.95    1.13
areaSI_EK:year_f1952     0.07    1.04
areaFB:year_f1953        0.15    1.06
areaFM:year_f1953        0.07    1.08
areaHO:year_f1953       -0.31    1.15
areaJM:year_f1953       -0.02    1.08
areaMU:year_f1953       -0.12    1.02
areaRA:year_f1953       -1.16    1.12
areaSI_EK:year_f1953    -0.30    1.04
areaFB:year_f1954        0.17    1.07
areaFM:year_f1954        0.08    1.11
areaHO:year_f1954       -0.16    1.17
areaJM:year_f1954       -0.10    1.09
areaMU:year_f1954       -0.40    1.05
areaRA:year_f1954       -1.36    1.11
areaSI_EK:year_f1954    -0.67    1.06
areaFB:year_f1955        0.70    1.08
areaFM:year_f1955        0.80    1.11
areaHO:year_f1955        0.26    1.17
areaJM:year_f1955        0.64    1.12
areaMU:year_f1955        0.10    1.05
areaRA:year_f1955       -1.51    1.12
areaSI_EK:year_f1955     0.07    1.08
areaFB:year_f1956        0.04    1.08
areaFM:year_f1956       -0.09    1.12
areaHO:year_f1956       -0.05    1.20
areaJM:year_f1956       -0.30    1.10
areaMU:year_f1956       -0.63    1.05
areaRA:year_f1956       -1.06    1.12
areaSI_EK:year_f1956    -0.94    1.07
areaFB:year_f1957        0.41    1.09
areaFM:year_f1957        0.41    1.12
areaHO:year_f1957        0.33    1.16
areaJM:year_f1957        0.30    1.11
areaMU:year_f1957       -0.26    1.06
areaRA:year_f1957       -0.91    1.12
areaSI_EK:year_f1957    -0.35    1.07
areaFB:year_f1958        0.45    1.15
areaFM:year_f1958        0.41    1.19
areaHO:year_f1958        0.50    1.20
areaJM:year_f1958        0.17    1.19
areaMU:year_f1958       -0.16    1.12
areaRA:year_f1958       -0.81    1.16
areaSI_EK:year_f1958    -0.27    1.16
areaFB:year_f1959        0.30    1.04
areaFM:year_f1959        0.24    1.08
areaHO:year_f1959       -0.15    1.12
areaJM:year_f1959        0.27    1.07
areaMU:year_f1959       -0.18    1.02
areaRA:year_f1959       -1.64    1.10
areaSI_EK:year_f1959    -0.17    1.04
areaFB:year_f1960        0.25    1.03
areaFM:year_f1960        0.22    1.08
areaHO:year_f1960       -0.13    1.12
areaJM:year_f1960        0.13    1.06
areaMU:year_f1960       -0.19    1.01
areaRA:year_f1960       -1.03    1.09
areaSI_EK:year_f1960    -0.28    1.04
areaFB:year_f1961        0.03    1.11
areaFM:year_f1961       -0.16    1.15
areaHO:year_f1961       -0.61    1.18
areaJM:year_f1961       -0.15    1.13
areaMU:year_f1961       -0.24    1.08
areaRA:year_f1961       -1.70    1.14
areaSI_EK:year_f1961    -0.49    1.09
areaFB:year_f1962        0.17    1.16
areaFM:year_f1962       -0.06    1.20
areaHO:year_f1962       -0.07    1.22
areaJM:year_f1962       -0.16    1.19
areaMU:year_f1962       -0.46    1.13
areaRA:year_f1962       -1.16    1.15
areaSI_EK:year_f1962    -0.95    1.16
areaFB:year_f1963        0.13    1.05
areaFM:year_f1963        0.01    1.10
areaHO:year_f1963       -0.24    1.15
areaJM:year_f1963       -0.20    1.09
areaMU:year_f1963       -0.30    1.01
areaRA:year_f1963       -1.07    1.13
areaSI_EK:year_f1963    -0.65    1.05
areaFB:year_f1964        0.06    1.10
areaFM:year_f1964       -0.17    1.14
areaHO:year_f1964       -0.05    1.18
areaJM:year_f1964       -0.27    1.11
areaMU:year_f1964       -0.42    1.07
areaRA:year_f1964       -0.92    1.15
areaSI_EK:year_f1964    -0.73    1.08
areaFB:year_f1965        0.22    1.10
areaFM:year_f1965        0.06    1.15
areaHO:year_f1965        0.22    1.18
areaJM:year_f1965        0.03    1.14
areaMU:year_f1965       -0.50    1.07
areaRA:year_f1965       -0.87    1.13
areaSI_EK:year_f1965    -0.85    1.10
areaFB:year_f1966        0.10    1.04
areaFM:year_f1966        0.05    1.08
areaHO:year_f1966       -0.37    1.14
areaJM:year_f1966       -0.04    1.06
areaMU:year_f1966       -0.11    1.01
areaRA:year_f1966       -1.24    1.10
areaSI_EK:year_f1966    -0.19    1.04
areaFB:year_f1967        0.26    1.08
areaFM:year_f1967        0.16    1.12
areaHO:year_f1967        0.03    1.16
areaJM:year_f1967        0.07    1.12
areaMU:year_f1967       -0.11    1.03
areaRA:year_f1967       -1.07    1.14
areaSI_EK:year_f1967    -0.23    1.06
areaFB:year_f1968        0.27    1.05
areaFM:year_f1968        0.19    1.10
areaHO:year_f1968        0.14    1.16
areaJM:year_f1968        0.14    1.09
areaMU:year_f1968       -0.21    1.02
areaRA:year_f1968       -1.19    1.13
areaSI_EK:year_f1968    -0.33    1.04
areaFB:year_f1969        0.40    1.03
areaFM:year_f1969        0.39    1.08
areaHO:year_f1969        0.02    1.13
areaJM:year_f1969        0.43    1.08
areaMU:year_f1969        0.06    1.00
areaRA:year_f1969       -1.02    1.10
areaSI_EK:year_f1969     0.24    1.04
areaFB:year_f1970        0.12    1.03
areaFM:year_f1970        0.04    1.08
areaHO:year_f1970       -0.09    1.13
areaJM:year_f1970       -0.02    1.06
areaMU:year_f1970       -0.22    1.00
areaRA:year_f1970       -0.85    1.11
areaSI_EK:year_f1970    -0.28    1.03
areaFB:year_f1971        0.23    1.06
areaFM:year_f1971        0.19    1.10
areaHO:year_f1971       -0.25    1.15
areaJM:year_f1971        0.30    1.09
areaMU:year_f1971       -0.03    1.03
areaRA:year_f1971       -1.25    1.09
areaSI_EK:year_f1971     0.06    1.05
areaFB:year_f1972        0.50    1.05
areaFM:year_f1972        0.52    1.11
areaHO:year_f1972       -0.32    1.13
areaJM:year_f1972        0.36    0.96
areaMU:year_f1972        0.39    1.01
areaRA:year_f1972       -1.01    1.11
areaSI_EK:year_f1972     0.41    1.04
areaFB:year_f1973        0.22    1.05
areaFM:year_f1973        0.18    1.08
areaHO:year_f1973       -0.01    1.13
areaJM:year_f1973       -2.75    0.98
areaMU:year_f1973       -0.22    1.05
areaRA:year_f1973       -0.94    1.10
areaSI_EK:year_f1973    -0.14    1.05
areaFB:year_f1974        0.34    1.06
areaFM:year_f1974        0.23    1.09
areaHO:year_f1974        0.15    1.15
areaJM:year_f1974       -3.73    0.93
areaMU:year_f1974       -0.06    1.02
areaRA:year_f1974       -1.25    1.11
areaSI_EK:year_f1974    -0.29    1.04
areaFB:year_f1975        0.51    1.05
areaFM:year_f1975       -0.65    0.99
areaHO:year_f1975        0.51    1.14
areaJM:year_f1975       -1.91    0.94
areaMU:year_f1975        0.02    1.03
areaRA:year_f1975       -0.60    1.10
areaSI_EK:year_f1975     0.11    1.05
areaFB:year_f1976       -0.45    0.97
areaFM:year_f1976        0.69    0.98
areaHO:year_f1976        0.53    1.14
areaJM:year_f1976        0.01    0.94
areaMU:year_f1976        0.05    1.02
areaRA:year_f1976       -0.43    1.11
areaSI_EK:year_f1976     0.36    1.04
areaFB:year_f1977       -2.00    0.94
areaFM:year_f1977       -0.77    0.96
areaHO:year_f1977       -0.07    1.17
areaJM:year_f1977       -1.53    0.93
areaMU:year_f1977        0.07    1.03
areaRA:year_f1977       -0.70    1.14
areaSI_EK:year_f1977     0.10    1.04
areaFB:year_f1978       -1.29    0.93
areaFM:year_f1978       -0.85    1.00
areaHO:year_f1978        0.00    1.14
areaJM:year_f1978       -1.13    0.94
areaMU:year_f1978        0.00    1.03
areaRA:year_f1978       -0.97    1.09
areaSI_EK:year_f1978     0.00    1.05
areaFB:year_f1979       -0.19    0.94
areaFM:year_f1979       -0.27    0.99
areaHO:year_f1979       -0.15    1.13
areaJM:year_f1979       -2.75    0.94
areaMU:year_f1979        0.13    1.00
areaRA:year_f1979       -0.56    1.11
areaSI_EK:year_f1979     0.25    1.03
areaFB:year_f1980       -1.01    0.98
areaFM:year_f1980        0.40    1.02
areaHO:year_f1980        0.12    1.16
areaJM:year_f1980       -0.26    0.94
areaMU:year_f1980        0.03    1.00
areaRA:year_f1980       -0.61    1.12
areaSI_EK:year_f1980     0.12    1.03
areaFB:year_f1981       -0.95    0.95
areaFM:year_f1981        0.22    1.01
areaHO:year_f1981       -0.35    1.15
areaJM:year_f1981       -0.69    0.93
areaMU:year_f1981        0.14    1.01
areaRA:year_f1981       -1.27    1.12
areaSI_EK:year_f1981     0.16    1.04
areaFB:year_f1982       -1.49    0.94
areaFM:year_f1982       -0.36    0.99
areaHO:year_f1982       -0.52    1.17
areaJM:year_f1982       -2.43    0.94
areaMU:year_f1982        0.11    1.02
areaRA:year_f1982       -0.66    1.14
areaSI_EK:year_f1982     0.25    1.05
areaFB:year_f1983       -1.37    0.92
areaFM:year_f1983        0.11    0.93
areaHO:year_f1983       -0.14    1.12
areaJM:year_f1983       -2.31    0.93
areaMU:year_f1983        0.01    1.01
areaRA:year_f1983       -0.95    1.09
areaSI_EK:year_f1983     0.10    1.03
areaFB:year_f1984       -0.05    0.96
areaFM:year_f1984        1.30    0.93
areaHO:year_f1984       -0.36    1.12
areaJM:year_f1984       -0.64    0.93
areaMU:year_f1984        0.02    1.01
areaRA:year_f1984       -1.27    1.10
areaSI_EK:year_f1984     0.03    1.05
areaFB:year_f1985        0.58    0.91
areaFM:year_f1985        3.93    0.93
areaHO:year_f1985        0.06    1.14
areaJM:year_f1985       -4.40    0.93
areaMU:year_f1985        0.11    1.01
areaRA:year_f1985       -0.44    1.13
areaSI_EK:year_f1985     0.29    1.03
areaFB:year_f1986       -1.36    0.93
areaFM:year_f1986        9.71    0.93
areaHO:year_f1986       -0.20    1.11
areaJM:year_f1986       -3.27    0.98
areaMU:year_f1986       -0.08    1.01
areaRA:year_f1986       -0.23    1.08
areaSI_EK:year_f1986    -0.17    1.04
areaFB:year_f1987       -1.06    0.93
areaFM:year_f1987        5.31    0.94
areaHO:year_f1987       -0.20    1.21
areaJM:year_f1987       -0.39    0.92
areaMU:year_f1987        0.06    1.04
areaRA:year_f1987       -0.05    1.18
areaSI_EK:year_f1987    -1.53    0.93
areaFB:year_f1988       -2.15    0.92
areaFM:year_f1988        6.85    0.94
areaHO:year_f1988        0.17    1.12
areaJM:year_f1988       -2.42    0.92
areaMU:year_f1988       -0.08    1.02
areaRA:year_f1988       -1.21    1.10
areaSI_EK:year_f1988    -3.87    0.95
areaFB:year_f1989       -3.13    0.91
areaFM:year_f1989        0.01    0.92
areaHO:year_f1989       -1.63    0.96
areaJM:year_f1989       -1.48    0.91
areaMU:year_f1989       -0.19    1.00
areaRA:year_f1989       -1.17    1.09
areaSI_EK:year_f1989    -1.41    0.89
areaFB:year_f1990       -2.13    0.91
areaFM:year_f1990        2.68    0.93
areaHO:year_f1990       -1.23    0.96
areaJM:year_f1990        1.26    0.92
areaMU:year_f1990       -0.26    1.02
areaRA:year_f1990       -1.05    1.09
areaSI_EK:year_f1990    -0.95    0.90
areaFB:year_f1991        0.57    0.80
areaFM:year_f1991        3.54    0.85
areaHO:year_f1991       -0.36    0.86
areaJM:year_f1991       -0.27    0.82
areaMU:year_f1991       -0.01    0.81
areaRA:year_f1991        0.38    0.86
areaSI_EK:year_f1991    -0.61    0.86
areaFB:year_f1992        0.88    0.78
areaFM:year_f1992        1.29    0.81
areaHO:year_f1992        0.02    0.85
areaJM:year_f1992        1.37    0.85
areaMU:year_f1992        1.00    0.83
areaRA:year_f1992        1.59    0.88
areaSI_EK:year_f1992    -1.08    0.78
areaFB:year_f1993        0.58    0.77
areaFM:year_f1993        0.74    0.81
areaHO:year_f1993       -0.81    0.84
areaJM:year_f1993        0.79    0.79
areaMU:year_f1993       -0.11    0.83
areaRA:year_f1993       -0.52    1.00
areaSI_EK:year_f1993    -0.45    0.78
areaFB:year_f1994        0.87    0.77
areaFM:year_f1994        2.49    0.80
areaHO:year_f1994        1.05    0.84
areaJM:year_f1994        1.37    0.79
areaMU:year_f1994       -1.39    0.84
areaRA:year_f1994        0.09    0.85
areaSI_EK:year_f1994     0.59    0.78
areaFB:year_f1995       -0.20    0.78
areaFM:year_f1995        0.86    0.82
areaHO:year_f1995       -1.57    0.85
areaJM:year_f1995        0.92    0.80
areaMU:year_f1995       -0.33    0.83
areaRA:year_f1995       -4.61    0.88
areaSI_EK:year_f1995    -3.42    0.79
areaFB:year_f1996        1.25    0.78
areaFM:year_f1996        2.84    0.81
areaHO:year_f1996        0.21    0.85
areaJM:year_f1996        0.69    0.83
areaMU:year_f1996       -0.92    0.77
areaRA:year_f1996        0.02    0.86
areaSI_EK:year_f1996    -1.42    0.78
areaFB:year_f1997       -1.57    0.79
areaFM:year_f1997       -1.82    0.82
areaHO:year_f1997       -2.83    0.91
areaJM:year_f1997       -0.73    0.84
areaMU:year_f1997       -3.51    0.77
areaRA:year_f1997       -2.05    0.88
areaSI_EK:year_f1997    -2.78    0.78
areaFB:year_f1998        1.02    0.76
areaFM:year_f1998        1.29    0.82
areaHO:year_f1998       -0.21    0.83
areaJM:year_f1998        1.06    0.78
areaMU:year_f1998       -0.92    0.74
areaRA:year_f1998       -0.03    0.82
areaSI_EK:year_f1998    -1.21    0.77
areaFB:year_f1999        1.29    0.76
areaFM:year_f1999        0.68    0.81
areaHO:year_f1999       -3.02    0.89
areaJM:year_f1999        1.10    0.78
areaMU:year_f1999       -1.39    0.74
areaRA:year_f1999       -1.62    0.82
areaSI_EK:year_f1999    -2.09    0.76
areaFB:year_f2000        1.22    0.76
areaFM:year_f2000        1.02    0.79
areaHO:year_f2000       -0.29    0.83
areaJM:year_f2000        1.72    0.77
areaMU:year_f2000       -0.40    0.74
areaRA:year_f2000       -0.75    0.81
areaSI_EK:year_f2000    -0.61    0.76
areaFB:year_f2001        0.72    0.76
areaFM:year_f2001       -0.09    0.79
areaHO:year_f2001       -1.33    0.83
areaJM:year_f2001        0.43    0.77
areaMU:year_f2001       -2.00    0.74
areaRA:year_f2001       -0.23    0.81
areaSI_EK:year_f2001    -2.79    0.76
areaFB:year_f2002        0.60    0.77
areaFM:year_f2002        1.56    0.80
areaHO:year_f2002       -1.78    0.83
areaJM:year_f2002        0.65    0.78
areaMU:year_f2002       -2.88    0.74
areaRA:year_f2002       -1.03    0.83
areaSI_EK:year_f2002    -1.85    0.76
areaFB:year_f2003        0.06    0.77
areaFM:year_f2003        0.72    0.81
areaHO:year_f2003       -0.73    0.84
areaJM:year_f2003       -0.19    0.78
areaMU:year_f2003       -2.72    0.75
areaRA:year_f2003       -1.33    0.83
areaSI_EK:year_f2003    -2.46    0.77
areaFB:year_f2004        0.91    0.76
areaFM:year_f2004        0.84    0.80
areaHO:year_f2004       -0.31    0.83
areaJM:year_f2004        0.95    0.77
areaMU:year_f2004       -1.48    0.74
areaRA:year_f2004       -1.12    0.82
areaSI_EK:year_f2004    -1.26    0.76
areaFB:year_f2005        0.63    0.76
areaFM:year_f2005        0.09    0.80
areaHO:year_f2005       -1.09    0.83
areaJM:year_f2005        0.33    0.78
areaMU:year_f2005       -2.06    0.74
areaRA:year_f2005       -0.78    0.82
areaSI_EK:year_f2005    -2.29    0.79
areaFB:year_f2006        3.35    0.80
areaFM:year_f2006        1.45    0.83
areaHO:year_f2006        1.80    0.86
areaJM:year_f2006        3.89    0.81
areaMU:year_f2006        0.47    0.77
areaRA:year_f2006        3.31    0.96
areaSI_EK:year_f2006     1.66    0.79
areaFB:year_f2007        1.72    0.79
areaFM:year_f2007        1.08    0.82
areaHO:year_f2007        0.36    0.85
areaJM:year_f2007        1.83    0.80
areaMU:year_f2007       -0.12    0.77
areaRA:year_f2007        0.06    0.84
areaSI_EK:year_f2007    -0.08    0.78
areaFB:year_f2008        1.14    0.79
areaFM:year_f2008        0.88    0.82
areaHO:year_f2008       -0.78    0.85
areaJM:year_f2008        1.16    0.80
areaMU:year_f2008       -0.37    0.77
areaRA:year_f2008       -1.05    0.85
areaSI_EK:year_f2008    -0.27    0.78
areaFB:year_f2009        0.79    0.79
areaFM:year_f2009       -0.05    0.82
areaHO:year_f2009       -0.79    0.86
areaJM:year_f2009        0.10    0.80
areaMU:year_f2009       -1.43    0.77
areaRA:year_f2009        0.15    0.85
areaSI_EK:year_f2009    -2.43    0.79
areaFB:year_f2010        1.27    0.90
areaFM:year_f2010        0.82    0.93
areaHO:year_f2010       -0.09    0.96
areaJM:year_f2010        1.52    0.91
areaMU:year_f2010       -0.58    0.88
areaRA:year_f2010        0.85    0.95
areaSI_EK:year_f2010     0.16    0.97
areaFB:year_f2011        2.09    0.90
areaFM:year_f2011        1.65    0.93
areaHO:year_f2011        1.04    0.96
areaJM:year_f2011        1.64    0.92
areaMU:year_f2011       -0.26    0.89
areaRA:year_f2011        1.34    0.95
areaSI_EK:year_f2011    -0.61    0.90
areaFB:year_f2012        1.21    0.88
areaFM:year_f2012        1.09    0.91
areaHO:year_f2012       -0.17    0.94
areaJM:year_f2012        0.82    0.89
areaMU:year_f2012       -0.72    0.86
areaRA:year_f2012        0.19    0.99
areaSI_EK:year_f2012    -1.09    0.88
areaFB:year_f2013        1.75    0.90
areaFM:year_f2013        1.57    0.93
areaHO:year_f2013        0.55    0.96
areaJM:year_f2013        1.68    0.91
areaMU:year_f2013       -0.13    0.88
areaRA:year_f2013       -0.15    1.00
areaSI_EK:year_f2013    -0.93    0.89
areaFB:year_f2014        1.63    0.89
areaFM:year_f2014        0.96    0.92
areaHO:year_f2014        0.13    0.95
areaJM:year_f2014        2.12    0.90
areaMU:year_f2014        0.20    0.87
areaRA:year_f2014        1.08    0.94
areaSI_EK:year_f2014    -0.18    0.88
areaFB:year_f2015        0.93    0.88
areaFM:year_f2015        0.49    0.91
areaHO:year_f2015       -0.58    0.96
areaJM:year_f2015        1.27    0.90
areaMU:year_f2015       -0.22    0.87
areaRA:year_f2015       -0.53    0.94
areaSI_EK:year_f2015    -0.65    0.88
areaFB:year_f2016        1.00    0.89
areaFM:year_f2016       -0.31    0.92
areaHO:year_f2016       -0.88    0.95
areaJM:year_f2016        0.82    0.90
areaMU:year_f2016       -0.97    0.87
areaRA:year_f2016        0.13    0.94
areaSI_EK:year_f2016    -1.15    0.88
areaFB:year_f2017        0.71    0.88
areaFM:year_f2017        0.22    0.91
areaHO:year_f2017       -0.83    0.94
areaJM:year_f2017        1.05    0.89
areaMU:year_f2017       -0.55    0.86
areaRA:year_f2017       -0.77    0.93
areaSI_EK:year_f2017    -0.81    0.87
areaFB:year_f2018        1.40    0.89
areaFM:year_f2018        1.54    0.92
areaHO:year_f2018       -0.87    0.95
areaJM:year_f2018        1.31    0.90
areaMU:year_f2018       -1.41    0.87
areaRA:year_f2018       -0.22    0.94
areaSI_EK:year_f2018    -1.57    0.89
areaFB:year_f2019       -0.04    0.96
areaFM:year_f2019        0.83    0.92
areaHO:year_f2019       -0.57    0.95
areaJM:year_f2019        1.28    0.90
areaMU:year_f2019       -0.78    0.87
areaRA:year_f2019       -0.32    0.94
areaSI_EK:year_f2019    -0.31    0.88
areaFB:year_f2020       -0.12    0.94
areaFM:year_f2020       -0.59    0.93
areaHO:year_f2020       -0.74    1.08
areaJM:year_f2020        1.03    0.93
areaMU:year_f2020        0.44    0.88
areaRA:year_f2020       -0.37    1.00
areaSI_EK:year_f2020    -1.67    1.05
areaFB:year_f2021        0.24    1.03
areaFM:year_f2021        0.66    0.98
areaHO:year_f2021       -0.92    1.05
areaJM:year_f2021        1.31    0.93
areaMU:year_f2021        1.63    0.90
areaRA:year_f2021       -2.02    1.00
areaSI_EK:year_f2021    -0.55    1.04
areaFB:year_f2022        0.48    1.74
areaFM:year_f2022        0.62    1.42
areaHO:year_f2022       -0.91    1.50
areaJM:year_f2022        0.43    1.43
areaMU:year_f2022        1.77    1.35
areaRA:year_f2022       -1.18    1.45
areaSI_EK:year_f2022     0.56    1.56

Smooth terms:
                     Std. Dev.
sds(yday):areaBS)         0.90
sds(yday):areaFB)         1.04
sds(yday):areaFM)         0.95
sds(yday):areaHO)         1.17
sds(yday):areaJM)         1.30
sds(yday):areaMU)         0.88
sds(yday):areaRA)         1.43
sds(yday):areaSI_EK)      0.75

Dispersion parameter: 1.70
ML criterion at convergence: 147129.577

See ?tidy.sdmTMB to extract these values as a data frame.

Predict

# Make a new data frame and predict!
nd <- data.frame(expand.grid(yday = seq(min(df$yday), max(df$yday), by = 1),
                             area = unique(df$area),
                             source = unique(df$source_f),
                             year = unique(df$year))) |>
  mutate(source_f = as.factor(source),
         year_f = as.factor(year)) |> # Left join in growth data column
  left_join(gdat, by = "area") |> 
  mutate(area = as.factor(area),
         growth_dat = ifelse(year >= min & year <= max, "Y", "N"))
mutate: new variable 'source_f' (factor) with 3 unique values and 0% NA
        new variable 'year_f' (factor) with 83 unique values and 0% NA
left_join: added 2 columns (min, max)
           > rows only in x         0
           > rows only in y  (      4)
           > matched rows     729,072
           >                 =========
           > rows total       729,072
mutate: converted 'area' from character to factor (0 new NA)
        new variable 'growth_dat' (character) with 2 unique values and 0% NA
# Predict
nd$pred <- predict(m, newdata = nd)$est

In order to have a lower temperature in before-nuclear times (without any data to inform that), we can use the nearby areas.. so FM informs BT prior to nuclear

nd_sub <- nd |> 
  mutate(keep = "N",
         keep = ifelse(area == "FM" & year <= 1980, "Y", keep), # use FM instead of BT
         keep = ifelse(area == "SI_EK" & year <= 1972, "Y", keep)) |> # use SI_EK instead of SI_HA
  filter(keep == "Y") |> # Now change the labels to BT and SI_EK...
  mutate(area = ifelse(area == "FM", "BT", "SI_HA"))
mutate: new variable 'keep' (character) with 2 unique values and 0% NA
filter: removed 647,820 rows (89%), 81,252 rows remaining
mutate: converted 'area' from factor to character (0 new NA)
# Bind rows and plot only the temperature series we will use for growth modelling
nd <- bind_rows(nd, nd_sub) |>
  select(-keep) |> 
  mutate(growth_dat = ifelse(area == "SI_HA" & year %in% c(1966, 1967), "Y", growth_dat)) # SI_EK and SI_HA do not have the same starting years, so we can't use allo columns from SI_EK
select: dropped one variable (keep)
mutate: changed 2,196 values (<1%) of 'growth_dat' (0 new NA)

Plot predictions

# Trimmed plot
nd |> 
  filter(growth_dat == "Y") |> 
  ggplot(aes(yday, y = year, fill = pred, color = pred)) +
  geom_raster() +
  facet_wrap(~area, ncol = 4) +
  scale_fill_viridis(option = "magma") +
  scale_color_viridis(option = "magma") +
  labs(x = "Day-of-the-year", y = "Year", color = "Predicted SST (°C)", fill = "Predicted SST (°C)") + 
  theme(legend.position = c(0.78, 0.14))
filter: removed 461,160 rows (57%), 349,164 rows remaining

ggsave(paste0(home, "/figures/supp/temp_pred_yday_area.pdf"), width = 17, height = 17, units = "cm")

# Full plot
nd |>
  ggplot(aes(yday, y = year, fill = pred, color = pred)) +
    geom_raster() +
    facet_wrap(~area, ncol = 4) +
    scale_fill_viridis(option = "magma") +
    scale_color_viridis(option = "magma") +
    labs(x = "Yearday", y = "Year", color = "Predicted SST (°C)", fill = "Predicted SST (°C)")

Detailed exploration of predictions

# Loop trough all areas, plot temperature as a function of yday, color by data source, facet by year

for(i in unique(nd$area)) {
  
  plotdat <- nd |> filter(area == i)
  
  print(
    ggplot(plotdat, aes(yday, pred, color = source)) + 
      scale_color_brewer(palette = "Dark2") + 
      facet_wrap(~year) + 
      geom_point(data = filter(d, area == i & year > min(plotdat$year)), size = 0.2,
                 aes(yday, temp, color = source)) + 
      geom_line(linewidth = 0.3) + 
      labs(title = paste("Area = ", i), color = "", linetype = "") + 
      guides(color = guide_legend(title.position = "top", title.hjust = 0.5)) + 
      theme_sleek(base_size = 8) +
      theme(legend.position = c(0.8, 0.08), 
            legend.direction = "horizontal",
            legend.spacing.y = unit(-0.3, "cm")) + 
      labs(x = "Day of the year", y = "Predicted SST (°C)")
  )
  
  ggsave(paste0(home, "/figures/supp/full_model/temp_pred_yday_area_full_", i, ".pdf" ), width = 17, height = 17, units = "cm")
  
}
filter: removed 719,190 rows (89%), 91,134 rows remaining
filter: removed 88,476 rows (96%), 3,449 rows remaining
filter: removed 719,190 rows (89%), 91,134 rows remaining
filter: removed 82,464 rows (90%), 9,461 rows remaining

filter: removed 719,190 rows (89%), 91,134 rows remaining
filter: removed 79,969 rows (87%), 11,956 rows remaining

filter: removed 719,190 rows (89%), 91,134 rows remaining
filter: removed 83,775 rows (91%), 8,150 rows remaining

filter: removed 719,190 rows (89%), 91,134 rows remaining
filter: removed 81,803 rows (89%), 10,122 rows remaining

filter: removed 719,190 rows (89%), 91,134 rows remaining
filter: removed 81,052 rows (88%), 10,873 rows remaining

filter: removed 719,190 rows (89%), 91,134 rows remaining
filter: removed 87,117 rows (95%), 4,808 rows remaining

filter: removed 719,190 rows (89%), 91,134 rows remaining
filter: removed 83,064 rows (90%), 8,861 rows remaining

filter: removed 765,306 rows (94%), 45,018 rows remaining
filter: removed 82,953 rows (90%), 8,972 rows remaining

filter: removed 774,090 rows (96%), 36,234 rows remaining
filter: removed 76,748 rows (83%), 15,177 rows remaining

Area-specific models

spec_preds <- list()
res_list <- list()

for(i in unique(d$area)) {
  
  dd <- d |> filter(area == i)

  if(unique(dd$area) %in% c("BS", "BT", "FB", "FM", "MU", "SI_EK")) { # RA, JM, HO, SI_HA remains
    
    mspec <- sdmTMB(temp ~ 0 + source_f + year_f + s(yday, bs = "cc"), 
                    data = dd,
                    family = student(df = 6),
                    spatial = "off",
                    spatiotemporal = "off",
                    knots = list(yday = c(0.5, 364.5)),
                    control = sdmTMBcontrol(newton_loops = 1)) 
  
  } else {
    
    mspec <- sdmTMB(temp ~ 0 + source_f + year_f + s(yday, bs = "cc"), 
                    data = dd,
                    family = student(df = 10),
                    spatial = "off",
                    spatiotemporal = "off",
                    knots = list(yday = c(0.5, 364.5)),
                    control = sdmTMBcontrol(newton_loops = 1)) 
    
  }
  
  sanity(mspec)

  # Residuals
  mcmc_res_msep <- residuals(mspec, type = "mle-mcmc",
                             mcmc_samples = sdmTMBextra::predict_mle_mcmc(mspec,
                                                                          mcmc_iter = 201,
                                                                          mcmc_warmup = 200))
  
  dd$res <- as.vector(mcmc_res_msep)
    
  # Store residuals
  res_list[[i]] <- dd
  
  print(ggplot(dd, aes(sample = mcmc_res_msep)) +
    stat_qq(size = 0.75, shape = 21, fill = NA) +
    stat_qq_line() +
    labs(y = "Sample Quantiles", x = "Theoretical Quantiles", title = paste("Area = ", i)) + 
    theme(aspect.ratio = 1))
  
  # Predict on new data
  nd_area <- data.frame(expand.grid(yday = seq(min(dd$yday), max(dd$yday), by = 1),
                                    area = unique(dd$area),
                                    source = unique(dd$source_f),
                                    year = unique(dd$year))) |>
    mutate(source_f = as.factor(source),
           year_f = as.factor(year)) |> 
    left_join(gdat, by = "area") |> 
    mutate(area = as.factor(area),
           growth_dat = ifelse(year >= min & year <= max, "Y", "N"))
  
  nd_area$pred <- predict(mspec, newdata = nd_area)$est

  # Save!
  spec_preds[[i]] <- nd_area
  
}
filter: removed 88,464 rows (96%), 3,461 rows remaining
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf

SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.000356 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 3.56 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 1.86447 seconds (Warm-up)
Chain 1:                0.002701 seconds (Sampling)
Chain 1:                1.86717 seconds (Total)
Chain 1: 
mutate: new variable 'source_f' (factor) with 3 unique values and 0% NA
        new variable 'year_f' (factor) with 83 unique values and 0% NA
left_join: added 2 columns (min, max)
           > rows only in x        0
           > rows only in y  (    11)
           > matched rows     83,664
           >                 ========
           > rows total       83,664
mutate: converted 'area' from character to factor (0 new NA)
        new variable 'growth_dat' (character) with 2 unique values and 0% NA
filter: removed 82,953 rows (90%), 8,972 rows remaining
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf

Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf


SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.000894 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 8.94 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 4.1908 seconds (Warm-up)
Chain 1:                0.003475 seconds (Sampling)
Chain 1:                4.19427 seconds (Total)
Chain 1: 
mutate: new variable 'source_f' (factor) with 3 unique values and 0% NA
        new variable 'year_f' (factor) with 42 unique values and 0% NA
left_join: added 2 columns (min, max)
           > rows only in x        0
           > rows only in y  (    11)
           > matched rows     46,116
           >                 ========
           > rows total       46,116
mutate: converted 'area' from character to factor (0 new NA)
        new variable 'growth_dat' (character) with 2 unique values and 0% NA
filter: removed 82,452 rows (90%), 9,473 rows remaining
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf

Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf


SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.000905 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 9.05 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 4.85745 seconds (Warm-up)
Chain 1:                0.006996 seconds (Sampling)
Chain 1:                4.86445 seconds (Total)
Chain 1: 
mutate: new variable 'source_f' (factor) with 3 unique values and 0% NA
        new variable 'year_f' (factor) with 83 unique values and 0% NA
left_join: added 2 columns (min, max)
           > rows only in x        0
           > rows only in y  (    11)
           > matched rows     90,885
           >                 ========
           > rows total       90,885
mutate: converted 'area' from character to factor (0 new NA)
        new variable 'growth_dat' (character) with 2 unique values and 0% NA
filter: removed 79,957 rows (87%), 11,968 rows remaining
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf

Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf


SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.001088 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 10.88 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 6.2292 seconds (Warm-up)
Chain 1:                0.008238 seconds (Sampling)
Chain 1:                6.23744 seconds (Total)
Chain 1: 
mutate: new variable 'source_f' (factor) with 3 unique values and 0% NA
        new variable 'year_f' (factor) with 83 unique values and 0% NA
left_join: added 2 columns (min, max)
           > rows only in x        0
           > rows only in y  (    11)
           > matched rows     91,134
           >                 ========
           > rows total       91,134
mutate: converted 'area' from character to factor (0 new NA)
        new variable 'growth_dat' (character) with 2 unique values and 0% NA
filter: removed 83,763 rows (91%), 8,162 rows remaining
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf

Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf


SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.000821 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 8.21 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 3.63301 seconds (Warm-up)
Chain 1:                0.012444 seconds (Sampling)
Chain 1:                3.64545 seconds (Total)
Chain 1: 
mutate: new variable 'source_f' (factor) with 3 unique values and 0% NA
        new variable 'year_f' (factor) with 83 unique values and 0% NA
left_join: added 2 columns (min, max)
           > rows only in x        0
           > rows only in y  (    11)
           > matched rows     90,885
           >                 ========
           > rows total       90,885
mutate: converted 'area' from character to factor (0 new NA)
        new variable 'growth_dat' (character) with 2 unique values and 0% NA
filter: removed 81,791 rows (89%), 10,134 rows remaining
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf

Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf


SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.000981 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 9.81 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 5.30214 seconds (Warm-up)
Chain 1:                0.007487 seconds (Sampling)
Chain 1:                5.30963 seconds (Total)
Chain 1: 
mutate: new variable 'source_f' (factor) with 3 unique values and 0% NA
        new variable 'year_f' (factor) with 83 unique values and 0% NA
left_join: added 2 columns (min, max)
           > rows only in x        0
           > rows only in y  (    11)
           > matched rows     84,660
           >                 ========
           > rows total       84,660
mutate: converted 'area' from character to factor (0 new NA)
        new variable 'growth_dat' (character) with 2 unique values and 0% NA
filter: removed 81,040 rows (88%), 10,885 rows remaining
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf

Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf


SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.00108 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 10.8 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 5.44634 seconds (Warm-up)
Chain 1:                0.00826 seconds (Sampling)
Chain 1:                5.4546 seconds (Total)
Chain 1: 
mutate: new variable 'source_f' (factor) with 3 unique values and 0% NA
        new variable 'year_f' (factor) with 83 unique values and 0% NA
left_join: added 2 columns (min, max)
           > rows only in x        0
           > rows only in y  (    11)
           > matched rows     91,134
           >                 ========
           > rows total       91,134
mutate: converted 'area' from character to factor (0 new NA)
        new variable 'growth_dat' (character) with 2 unique values and 0% NA
filter: removed 87,105 rows (95%), 4,820 rows remaining
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf

Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf


SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.000508 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 5.08 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 2.53892 seconds (Warm-up)
Chain 1:                0.003867 seconds (Sampling)
Chain 1:                2.54279 seconds (Total)
Chain 1: 
mutate: new variable 'source_f' (factor) with 3 unique values and 0% NA
        new variable 'year_f' (factor) with 83 unique values and 0% NA
left_join: added 2 columns (min, max)
           > rows only in x        0
           > rows only in y  (    11)
           > matched rows     91,134
           >                 ========
           > rows total       91,134
mutate: converted 'area' from character to factor (0 new NA)
        new variable 'growth_dat' (character) with 2 unique values and 0% NA
filter: removed 83,052 rows (90%), 8,873 rows remaining
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf

Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf


SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.000919 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 9.19 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 4.79943 seconds (Warm-up)
Chain 1:                0.006822 seconds (Sampling)
Chain 1:                4.80625 seconds (Total)
Chain 1: 
mutate: new variable 'source_f' (factor) with 3 unique values and 0% NA
        new variable 'year_f' (factor) with 83 unique values and 0% NA
left_join: added 2 columns (min, max)
           > rows only in x        0
           > rows only in y  (    11)
           > matched rows     91,134
           >                 ========
           > rows total       91,134
mutate: converted 'area' from character to factor (0 new NA)
        new variable 'growth_dat' (character) with 2 unique values and 0% NA
filter: removed 76,748 rows (83%), 15,177 rows remaining
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf

Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf


SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.001461 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 14.61 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 7.90052 seconds (Warm-up)
Chain 1:                0.005614 seconds (Sampling)
Chain 1:                7.90614 seconds (Total)
Chain 1: 
mutate: new variable 'source_f' (factor) with 3 unique values and 0% NA
        new variable 'year_f' (factor) with 50 unique values and 0% NA
left_join: added 2 columns (min, max)
           > rows only in x        0
           > rows only in y  (    11)
           > matched rows     54,900
           >                 ========
           > rows total       54,900
mutate: converted 'area' from character to factor (0 new NA)
        new variable 'growth_dat' (character) with 2 unique values and 0% NA

nd_area <- dplyr::bind_rows(spec_preds)
area_res <- dplyr::bind_rows(res_list)

# Plot residuals
dfs <- data.frame(area = c("BS", "BT", "FB", "FM", "MU", "SI_EK", "RA", "JM", "HO", "SI_HA"),
                  df = c(rep("\u03BD = 6", 7), rep("\u03BD = 10", 3)))
  
area_res |> 
  ggplot(aes(sample = res)) +
  stat_qq(size = 0.75, shape = 21, fill = NA) +
  facet_wrap(~area) +
  stat_qq_line() +
  labs(y = "Sample Quantiles", x = "Theoretical Quantiles") +
  geom_text(data = dfs, aes(label = df, hjust = -0.1, vjust = 1.25), inherit.aes = FALSE,
              x = -Inf, y = Inf, size = 2.6, color = "grey20") +
  theme(aspect.ratio = 1)

ggsave(paste0(home, "/figures/supp/qq_temp_area.pdf"), width = 17, height = 17, units = "cm", device = cairo_pdf)

# In order to have a lower temperature in before-nuclear times (without any data to inform that), we can use the nearby areas.. so FM informs BT prior to nuclear
nd_area_sub <- nd_area |> 
  mutate(keep = "N",
         keep = ifelse(area == "FM" & year <= 1980, "Y", keep), # use FM instead of BT
         keep = ifelse(area == "SI_EK" & year <= 1972, "Y", keep)) |> # use SI_EK instead of SI_HA
  filter(keep == "Y") |> # Now change the labels to BT and SI_EK...
  mutate(area = ifelse(area == "FM", "BT", "SI_HA"))
mutate: new variable 'keep' (character) with 2 unique values and 0% NA
filter: removed 734,394 rows (90%), 81,252 rows remaining
mutate: converted 'area' from factor to character (0 new NA)
# Bind rows and plot only the temperature series we will use for growth modelling
nd_area <- bind_rows(nd_area, nd_area_sub) |>
  select(-keep) |> 
  mutate(growth_dat = ifelse(area == "SI_HA" & year %in% c(1966, 1967), "Y", growth_dat)) # SI_EK and SI_HA do not have the same starting years, so we can't use allo columns from SI_EK
select: dropped one variable (keep)
mutate: changed 2,196 values (<1%) of 'growth_dat' (0 new NA)
nd_area |> 
  filter(area %in% c("FM", "BT", "SI_EK", "SI_HA")) |> 
  filter(year <= 1980 & year >= 1966) |> 
  group_by(area, year) |> 
  summarise(mean_temp = mean(pred)) |> 
  ungroup() |> 
  pivot_wider(names_from = area, values_from = mean_temp)
filter: removed 532,362 rows (59%), 364,536 rows remaining
filter: removed 298,656 rows (82%), 65,880 rows remaining
group_by: 2 grouping variables (area, year)
summarise: now 60 rows and 3 columns, one group variable remaining (area)
ungroup: no grouping variables
pivot_wider: reorganized (area, mean_temp) into (BT, FM, SI_EK, SI_HA) [was 60x3, now 15x5]
# A tibble: 15 × 5
    year    BT    FM SI_EK SI_HA
   <dbl> <dbl> <dbl> <dbl> <dbl>
 1  1966  8.15  8.15  7.97  7.97
 2  1967  8.98  8.98  8.76  8.76
 3  1968  8.69  8.69  8.37  8.37
 4  1969  8.55  8.55  8.53  8.53
 5  1970  8.30  8.30  8.04  8.04
 6  1971  8.33  8.33  8.33  8.33
 7  1972  8.83  8.83  8.90  8.90
 8  1973  9.19  9.19  8.99 11.1 
 9  1974  9.05  9.05  8.66  8.46
10  1975  7.43  7.43  9.34 14.6 
11  1976  7.09  7.09  8.00 13.7 
12  1977  5.61  5.61  7.80 13.0 
13  1978  5.59  5.59  7.57 13.1 
14  1979  5.86  5.86  7.42 11.4 
15  1980  6.96  6.96  7.73 13.4 

Plot detailed exploration of predictions

# Loop trough all areas, plot temperature as a function of yday, color by data source, facet by year

for(i in unique(nd_area$area)) {
  
  plotdat <- nd_area |> filter(area == i)
  
  print(
    ggplot(plotdat, aes(yday, pred, color = source)) + 
      scale_color_brewer(palette = "Dark2") + 
      facet_wrap(~year) + 
      geom_point(data = filter(d, area == i & year > min(plotdat$year)), size = 0.2,
                 aes(yday, temp, color = source)) + 
      geom_line(linewidth = 0.3) + 
      labs(title = paste("Area = ", i), color = "", linetype = "") + 
      guides(color = guide_legend(title.position = "top", title.hjust = 0.5)) + 
      theme_sleek(base_size = 8) +
      theme(legend.position = "bottom", 
            legend.direction = "horizontal",
            legend.spacing.y = unit(-0.3, "cm")) + 
      labs(x = "Day of the year", y = "Predicted SST (°C)")
  )
  
  ggsave(paste0(home, "/figures/supp/temp_pred_yday_area_", i, ".pdf" ), width = 17, height = 17, units = "cm")
  
}
filter: removed 813,234 rows (91%), 83,664 rows remaining
filter: removed 88,476 rows (96%), 3,449 rows remaining
filter: removed 805,764 rows (90%), 91,134 rows remaining
filter: removed 82,953 rows (90%), 8,972 rows remaining

filter: removed 806,013 rows (90%), 90,885 rows remaining
filter: removed 82,464 rows (90%), 9,461 rows remaining

filter: removed 805,764 rows (90%), 91,134 rows remaining
filter: removed 79,969 rows (87%), 11,956 rows remaining

filter: removed 806,013 rows (90%), 90,885 rows remaining
filter: removed 83,775 rows (91%), 8,150 rows remaining

filter: removed 812,238 rows (91%), 84,660 rows remaining
filter: removed 81,803 rows (89%), 10,122 rows remaining

filter: removed 805,764 rows (90%), 91,134 rows remaining
filter: removed 81,052 rows (88%), 10,873 rows remaining

filter: removed 805,764 rows (90%), 91,134 rows remaining
filter: removed 87,117 rows (95%), 4,808 rows remaining

filter: removed 805,764 rows (90%), 91,134 rows remaining
filter: removed 83,064 rows (90%), 8,861 rows remaining

filter: removed 805,764 rows (90%), 91,134 rows remaining
filter: removed 76,748 rows (83%), 15,177 rows remaining

# Same but trimmed
for(i in unique(nd_area$area)) {
  
  plotdat <- nd_area |> filter(area == i) |> filter(growth_dat == "Y")
  
  print(
    ggplot(plotdat, aes(yday, pred, color = source)) + 
      scale_color_brewer(palette = "Dark2") + 
      facet_wrap(~year) + 
      geom_point(data = filter(d, area == i & year > min(plotdat$year)), size = 0.2,
                 aes(yday, temp, color = source)) + 
      geom_line(linewidth = 0.3) + 
      labs(title = paste("Area = ", i), color = "", linetype = "") + 
      guides(color = guide_legend(title.position = "top", title.hjust = 0.5)) + 
      theme_sleek(base_size = 8) +
      theme(#legend.position = c(0.8, 0.08), 
            legend.position = "bottom", 
            legend.direction = "horizontal",
            legend.spacing.y = unit(-0.3, "cm")) + 
      labs(x = "Day of the year", y = "Predicted SST (°C)")
  )
  
  ggsave(paste0(home, "/figures/supp/temp_pred_yday_area_trimmed_", i, ".pdf" ), width = 17, height = 17, units = "cm")
  
}
filter: removed 813,234 rows (91%), 83,664 rows remaining
filter: removed 65,520 rows (78%), 18,144 rows remaining
filter: removed 89,016 rows (97%), 2,909 rows remaining
filter: removed 805,764 rows (90%), 91,134 rows remaining
filter: removed 49,410 rows (54%), 41,724 rows remaining
filter: removed 82,953 rows (90%), 8,972 rows remaining

filter: removed 806,013 rows (90%), 90,885 rows remaining
filter: removed 38,325 rows (42%), 52,560 rows remaining
filter: removed 82,812 rows (90%), 9,113 rows remaining

filter: removed 805,764 rows (90%), 91,134 rows remaining
filter: removed 49,410 rows (54%), 41,724 rows remaining
filter: removed 80,245 rows (87%), 11,680 rows remaining

filter: removed 806,013 rows (90%), 90,885 rows remaining
filter: removed 52,560 rows (58%), 38,325 rows remaining
filter: removed 84,279 rows (92%), 7,646 rows remaining

filter: removed 812,238 rows (91%), 84,660 rows remaining
filter: removed 19,380 rows (23%), 65,280 rows remaining
filter: removed 81,959 rows (89%), 9,966 rows remaining

filter: removed 805,764 rows (90%), 91,134 rows remaining
filter: removed 69,174 rows (76%), 21,960 rows remaining
filter: removed 81,544 rows (89%), 10,381 rows remaining

filter: removed 805,764 rows (90%), 91,134 rows remaining
filter: removed 66,978 rows (73%), 24,156 rows remaining
filter: removed 87,657 rows (95%), 4,268 rows remaining

filter: removed 805,764 rows (90%), 91,134 rows remaining
filter: removed 38,430 rows (42%), 52,704 rows remaining
filter: removed 83,400 rows (91%), 8,525 rows remaining

filter: removed 805,764 rows (90%), 91,134 rows remaining
filter: removed 37,332 rows (41%), 53,802 rows remaining
filter: removed 76,748 rows (83%), 15,177 rows remaining

Compare predictions from full and area-specific models

long_preds <- bind_rows(nd_area |> filter(growth_dat == "Y") |> mutate(model = "area-specific"),
                        nd |> filter(growth_dat == "Y") |> mutate(model = "full model"))
filter: removed 486,519 rows (54%), 410,379 rows remaining
mutate: new variable 'model' (character) with one unique value and 0% NA
filter: removed 461,160 rows (57%), 349,164 rows remaining
mutate: new variable 'model' (character) with one unique value and 0% NA
for(i in unique(long_preds$area)) {
  
  plotdat <- long_preds |> filter(area == i)
  
  print(
    ggplot(plotdat, aes(yday, pred, color = source, linetype = model)) + 
      scale_color_brewer(palette = "Dark2") + 
      facet_wrap(~year) + 
      geom_point(data = filter(d, area == i & growth_dat == "Y"), size = 0.2,
                 aes(yday, temp, color = source), inherit.aes = FALSE) + 
      geom_line(linewidth = 0.3) + 
      guides(color = guide_legend(title.position = "top", title.hjust = 0.5)) + 
      theme_sleek(base_size = 8) +
      theme(legend.position = "bottom", 
            legend.direction = "horizontal",
            legend.spacing.y = unit(-0.3, "cm")) + 
      labs(x = "Day of the year", y = "Predicted SST (°C)",
           title = paste("Area = ", i), color = "", linetype = "")
  )
  
  ggsave(paste0(home, "/figures/supp/full_model/temp_pred_comp_area_model_", i, ".pdf" ), width = 17, height = 17, units = "cm")
  
}
filter: removed 721,635 rows (95%), 37,908 rows remaining
filter: removed 90,033 rows (98%), 1,892 rows remaining
filter: removed 698,055 rows (92%), 61,488 rows remaining
filter: removed 85,553 rows (93%), 6,372 rows remaining

filter: removed 654,279 rows (86%), 105,264 rows remaining
filter: removed 83,349 rows (91%), 8,576 rows remaining

filter: removed 676,095 rows (89%), 83,448 rows remaining
filter: removed 85,630 rows (93%), 6,295 rows remaining

filter: removed 682,788 rows (90%), 76,755 rows remaining
filter: removed 85,073 rows (93%), 6,852 rows remaining

filter: removed 623,991 rows (82%), 135,552 rows remaining
filter: removed 82,998 rows (90%), 8,927 rows remaining

filter: removed 715,623 rows (94%), 43,920 rows remaining
filter: removed 89,965 rows (98%), 1,960 rows remaining

filter: removed 711,231 rows (94%), 48,312 rows remaining
filter: removed 89,689 rows (98%), 2,236 rows remaining

filter: removed 654,135 rows (86%), 105,408 rows remaining
filter: removed 84,947 rows (92%), 6,978 rows remaining

filter: removed 698,055 rows (92%), 61,488 rows remaining
filter: removed 78,821 rows (86%), 13,104 rows remaining

# Join wide data to plot differences easily
wide_pred <- left_join(nd_area |> 
                         dplyr::select(yday, area, year, pred, source) |> 
                         pivot_wider(names_from = source, values_from = pred) |> 
                         rename(logger_as = logger,
                                errs_as = errs,
                                fishing_as = fishing),
                       
                       nd |> 
                         dplyr::select(yday, area, year, pred, source) |> 
                         pivot_wider(names_from = source, values_from = pred) |> 
                         rename(logger_full = logger,
                                errs_full = errs,
                                fishing_full = fishing),
                       
                       by = c("yday", "area", "year"))
pivot_wider: reorganized (pred, source) into (logger, errs, fishing) [was 896898x5, now 298966x6]
rename: renamed 3 variables (logger_as, errs_as, fishing_as)
pivot_wider: reorganized (pred, source) into (logger, errs, fishing) [was 810324x5, now 270108x6]
rename: renamed 3 variables (logger_full, errs_full, fishing_full)
left_join: added 3 columns (logger_full, errs_full, fishing_full)
           > rows only in x    33,672
           > rows only in y  (  4,814)
           > matched rows     265,294
           >                 =========
           > rows total       298,966
diff_pred <- wide_pred |> 
  mutate(logger_diff = logger_as - logger_full,
         errs_diff = errs_as - errs_full,
         fishing_diff = fishing_as - fishing_full) |> 
  select(yday, area, year, logger_diff, errs_diff, fishing_diff) |> 
  pivot_longer(c("logger_diff", "errs_diff", "fishing_diff"))
mutate: new variable 'logger_diff' (double) with 238,211 unique values and 11% NA
        new variable 'errs_diff' (double) with 238,211 unique values and 11% NA
        new variable 'fishing_diff' (double) with 238,211 unique values and 11% NA
select: dropped 6 variables (logger_as, errs_as, fishing_as, logger_full, errs_full, …)
pivot_longer: reorganized (logger_diff, errs_diff, fishing_diff) into (name, value) [was 298966x6, now 896898x5]
for(i in unique(diff_pred$area)) {
  
  plotdat <- diff_pred |> filter(area == i)
  
  print(ggplot(plotdat, aes(yday, value, color = name)) + 
      facet_wrap(~year) +  
      scale_color_brewer(palette = "Dark2") + 
      geom_hline(yintercept = 0, linetype = 2, alpha = 0.5) + 
      geom_line() +
      labs(x = "Day of the year", y = "Predicted SST (°C)",
           title = paste("Area = ", i), color = "", linetype = ""))
}
filter: removed 813,234 rows (91%), 83,664 rows remaining
filter: removed 805,764 rows (90%), 91,134 rows remaining

Warning: Removed 46116 rows containing missing values (`geom_line()`).
filter: removed 806,013 rows (90%), 90,885 rows remaining

filter: removed 805,764 rows (90%), 91,134 rows remaining

filter: removed 806,013 rows (90%), 90,885 rows remaining

filter: removed 812,238 rows (91%), 84,660 rows remaining

filter: removed 805,764 rows (90%), 91,134 rows remaining

filter: removed 805,764 rows (90%), 91,134 rows remaining

filter: removed 805,764 rows (90%), 91,134 rows remaining

filter: removed 805,764 rows (90%), 91,134 rows remaining

Warning: Removed 54900 rows containing missing values (`geom_line()`).

Plot summarized data and predictions

# Summarise data
dsum <- d |> 
  group_by(year, area, source) |> 
  summarise(temp = mean(temp)) |> 
  mutate(type = "data")
group_by: 3 grouping variables (year, area, source)
summarise: now 1,376 rows and 4 columns, 2 group variables remaining (year, area)
mutate (grouped): new variable 'type' (character) with one unique value and 0% NA
# Summarise predictions from full model and area-specific model
preds <- nd |> 
  filter(growth_dat == "Y" & source == "logger") |> 
  group_by(area, year) |> 
  summarise(temp = mean(pred)) |> 
  mutate(model = "full model")
filter: removed 693,936 rows (86%), 116,388 rows remaining
group_by: 2 grouping variables (area, year)
summarise: now 318 rows and 3 columns, one group variable remaining (area)
mutate (grouped): new variable 'model' (character) with one unique value and 0% NA
preds_area <- nd_area |> 
  filter(growth_dat == "Y" & source == "logger") |> 
  group_by(area, year) |> 
  summarise(temp = mean(pred)) |> 
  mutate(model = "area model")
filter: removed 760,105 rows (85%), 136,793 rows remaining
group_by: 2 grouping variables (area, year)
summarise: now 380 rows and 3 columns, one group variable remaining (area)
mutate (grouped): new variable 'model' (character) with one unique value and 0% NA
preds_comb <- bind_rows(preds, preds_area)

ggplot(preds_comb, aes(year, temp, color = source, linetype = model)) + 
  geom_point(data = dsum, aes(year, temp, color = source), size = 0.75, alpha = 0.75, inherit.aes = FALSE) + 
  scale_color_brewer(palette = "Accent") +
  geom_line(linewidth = 0.5, color = "grey20") + 
  facet_wrap(~area, ncol = 3) +
  theme(legend.position = "bottom")

Make final plot using the area-specific model

# Overall t_opt from 02-fit-vbge.qmd! Update when I got final model!
topt_lwr <- 6.402944
topt <- 9.468607
topt_upr <- 14.206451

# Add latitude
area <- c("BS", "BT", "FB", "FM", "HO", "JM", "MU", "RA", "SI_EK", "SI_HA")
nareas <- length(area)
lat <- c(60, 60.4, 60.3, 60.5, 63.7, 58, 59, 65.9, 57.3, 57.4)
lon <- c(21.5, 18.1, 19.5, 18, 20.9, 16.8, 18.1, 22.3, 16.6, 16.7)
area_attr <- data.frame(cbind(area = area, lat = lat, lon = lon)) |>
  mutate_at(c("lat","lon"), as.numeric) |> 
  arrange(desc(lat))
mutate_at: converted 'lat' from character to double (0 new NA)
           converted 'lon' from character to double (0 new NA)
ggplot(preds_area, aes(year, temp, color = temp)) + 
  facet_wrap(~factor(area, levels = area_attr$area), ncol = 4) + 
  geom_hline(yintercept = topt_lwr, linewidth = 0.3, linetype = 2, color = "tomato3", alpha = 0.3) +
  geom_hline(yintercept = topt, linewidth = 0.3, linetype = 1, color = "tomato3", alpha = 0.3) +
  geom_hline(yintercept = topt_upr, linewidth = 0.3, linetype = 2, color = "tomato3", alpha = 0.3) +
  geom_line(linewidth = 0.6) +
  labs(x = "Year", y = "Model-predicted annual average temperature") + 
  scale_color_viridis(option = "magma", name = "Area") +
  guides(color = "none") 

ggsave(paste0(home, "/figures/annual_average_temperature.pdf"), width = 17, height = 12, units = "cm")

# Check year range
preds_area |> 
  group_by(area) |> 
  summarise(min = min(year),
            max = max(year)) |> 
  arrange(min)
group_by: one grouping variable (area)
summarise: now 10 rows and 3 columns, ungrouped
# A tibble: 10 × 3
   area    min   max
   <chr> <dbl> <dbl>
 1 JM     1953  2016
 2 BT     1963  2000
 3 FM     1963  2000
 4 SI_HA  1966  2014
 5 SI_EK  1968  2015
 6 FB     1969  2016
 7 MU     1981  2000
 8 HO     1982  2016
 9 BS     1985  2002
10 RA     1985  2006
# Check temperature range
preds_area |> 
  group_by(area) |> 
  summarise(min = min(temp),
            max = max(temp)) |> 
  arrange(min)
group_by: one grouping variable (area)
summarise: now 10 rows and 3 columns, ungrouped
# A tibble: 10 × 3
   area    min   max
   <chr> <dbl> <dbl>
 1 RA     3.06  9.56
 2 JM     3.37 11.8 
 3 HO     3.99  8.07
 4 FB     5.04 10.6 
 5 MU     5.16  8.57
 6 SI_EK  5.49 10.4 
 7 BT     5.87 16.6 
 8 FM     5.87 15.0 
 9 BS     6.22 11.8 
10 SI_HA  7.76 16.9 
# Save prediction df
write_csv(preds_area, paste0(home, "/output/gam_predicted_temps.csv"))
knitr::knit_exit()

Code below is not evaluated!

Growing season? This might be different for different areas…

# Find day of the year where temperature exceeds 10C by area across years
# TODO: Also for area-specific predictions?
gs_area <- nd |> 
  group_by(area, yday) |> 
  summarise(mean_pred = mean(pred)) |>
  ungroup() |> 
  filter(mean_pred > 10) |> 
  group_by(area) |> 
  summarise(gs_min = min(yday),
            gs_max = max(yday))
group_by: 2 grouping variables (area, yday)
summarise: now 3,660 rows and 3 columns, one group variable remaining (area)
ungroup: no grouping variables
filter: removed 2,287 rows (62%), 1,373 rows remaining
group_by: one grouping variable (area)
summarise: now 10 rows and 3 columns, ungrouped
nd <- left_join(nd, gs_area, by = "area")
left_join: added 2 columns (gs_min, gs_max)
           > rows only in x         0
           > rows only in y  (      0)
           > matched rows     810,324
           >                 =========
           > rows total       810,324
gs_area$mean_pred <- 10

# Plot!
nd |> 
  #filter(growth_dat == "Y") |> 
  group_by(area, yday) |> 
  summarise(mean_pred = mean(pred)) |> 
  ggplot(aes(yday, mean_pred)) +
  geom_line() +
  labs(x = "Yearday", y = "Year", color = "Predicted SST (°C)", fill = "Predicted SST (°C)") + 
  facet_wrap(~area) +
  geom_point(data = gs_area, aes(gs_min, mean_pred), inherit.aes = FALSE, color = "tomato2") + 
  geom_point(data = gs_area, aes(gs_max, mean_pred), inherit.aes = FALSE, color = "tomato2") +
  geom_hline(yintercept = 10, linetype = 2)
group_by: 2 grouping variables (area, yday)
summarise: now 3,660 rows and 3 columns, one group variable remaining (area)

Now see if there is a systematic pattern in the difference between predicted and observed logger data, which could indicate that the source effect isn’t global but area-specific.

dlog <- d |> 
  filter(source == "logger") |> 
  mutate(type = "data",
         id = paste(area, year, yday, sep = "_")) |> 
  select(id, temp) |> 
  group_by(id) |> 
  summarise(obs = mean(temp)) # sometimes we have more than 1 observation per id
filter: removed 30,164 rows (33%), 61,761 rows remaining
mutate: changed 61,761 values (100%) of 'id' (61761 fewer NA)
        new variable 'type' (character) with one unique value and 0% NA
select: dropped 18 variables (year, area, yday, month, date, …)
group_by: one grouping variable (id)
summarise: now 58,744 rows and 2 columns, ungrouped
# dlog |> 
#   group_by(id) |> 
#   summarise(n = n()) |> 
#   distinct(n)

preds_log <- nd |> 
  filter(growth_dat == "Y" & source == "logger") |> 
  mutate(type = "model",
         id = paste(area, year, yday, sep = "_")) |> 
  filter(id %in% unique(dlog$id)) |> 
  ungroup() |> 
  left_join(dlog, by = "id")
filter: removed 693,936 rows (86%), 116,388 rows remaining
mutate: new variable 'type' (character) with one unique value and 0% NA
        new variable 'id' (character) with 116,388 unique values and 0% NA
filter: removed 89,556 rows (77%), 26,832 rows remaining
ungroup: no grouping variables
left_join: added one column (obs)
           > rows only in x        0
           > rows only in y  (31,912)
           > matched rows     26,832
           >                 ========
           > rows total       26,832
p1 <- preds_log |> 
  mutate(resid = pred - obs) |> 
  ggplot(aes(as.factor(area), resid, group = as.factor(area))) +
  #geom_jitter(alpha = 0.05, color = "grey20", height = 0, width = 0.2) + 
  geom_violin(fill = "grey70", color = NA) +
  geom_boxplot(width = 0.2, outlier.colour = NA, outlier.color = NA, outlier.fill = NA) +
  guides(color = "none") +
  geom_hline(yintercept = 0, linetype = 2, color = "tomato3", linewidth = 0.75) + 
  labs(x = "Area", y = "Manual residuals")
mutate: new variable 'resid' (double) with 26,832 unique values and 0% NA
p1

p1 + facet_wrap(~year) 

Some extra plots for BT especially, because it seems that the offset for logger isn’t that big in years without any logger data.. First predict BT only but with all data sources. Note this code is not run anymore and is instead expanded above to work on all areas

# Predict
ndbt <- data.frame(expand.grid(yday = seq(min(d$yday), max(d$yday), by = 1),
                               year = seq(1980, #filter(gdat, area == "BT")$min, 
                                          2004), #filter(gdat, area == "BT")$min),
                               source = unique(d$source))) |>
  mutate(area = "BT") |> 
  mutate(area = as.factor(area),
         source_f = as.factor(source),
         year_f = as.factor(year)) 

ndbt$pred <- predict(m, newdata = ndbt)$est

# Plot
ggplot(ndbt, aes(yday, pred, color = source)) + 
  scale_color_brewer(palette = "Accent") + 
  facet_wrap(~year) + 
  geom_point(data = filter(d, area == "BT" & year >= min(ndbt$year) & year <= max(ndbt$year)), size = 0.2,
             aes(yday, temp)) + 
  geom_line() + 
  labs(title = "Area = BT", color = "") + 
  guides(color = guide_legend(title.position = "top", title.hjust = 0.5)) + 
  theme_sleek(base_size = 8) +
  theme(legend.position = "bottom", 
        legend.direction = "horizontal",
        legend.spacing.y = unit(-0.3, "cm")) + 
  labs(x = "Day of the year", y = "Predicted SST (°C)")

ggsave(paste0(home, "/figures/supp/BT_test_plots.pdf" ), width = 17, height = 17, units = "cm")

# Right, so there is an offset... but it's not big. And the cold temperatures in
# certain years in the 1980's is not due to the offset not working, just that the year
# effect in that year is only informed by cold temperatures and there's no "memory"... 
# Though it seems also that there could be a bigger offset.. perhaps this indicates
# the source should indeed vary by area? Maybe, we if we fit models by area separately instead? Then we don't
# need the area interaction, and we can instead use a random walk or AR1 structure one
# the year effect?

dbt <- d |> filter(area == "BT")

mbt <- sdmTMB(temp ~ 0 + source_f + year_f + s(yday, bs = "cc"), 
              data = dbt,
              family = student(df = 5),
              spatial = "off",
              spatiotemporal = "off",
              knots = list(yday = c(0.5, 364.5)),
              control = sdmTMBcontrol(newton_loops = 1))

sanity(mbt)
summary(mbt)

# Predict
ndbt2 <- data.frame(expand.grid(yday = seq(min(d$yday), max(d$yday), by = 1),
                                year = seq(1980, #filter(gdat, area == "BT")$min, 
                                           2004), #filter(gdat, area == "BT")$min),
                                source = unique(d$source))) |>
  mutate(source_f = as.factor(source),
         year_f = as.factor(year)) 

ndbt2$pred <- predict(mbt, newdata = ndbt2)$est

# Plot without data
ggplot(filter(d, area == "BT" & year >= min(ndbt$year) & year <= max(ndbt$year)),
       aes(yday, temp, color = source)) + 
  geom_line(data = ndbt, aes(yday, pred, color = source, linetype = "main model"), alpha = 0.2) +
  geom_line(data = ndbt2, aes(yday, pred, color = source, linetype = "BT specific model"), alpha = 0.2) + 
  scale_color_brewer(palette = "Dark2") + 
  facet_wrap(~year) + 
  labs(title = "Area = BT", color = "", linetype = "") + 
  guides(color = guide_legend(title.position = "top", title.hjust = 0.5)) + 
  theme_sleek(base_size = 8) +
  theme(legend.position = "bottom", 
        legend.direction = "horizontal",
        legend.spacing.y = unit(-0.3, "cm")) + 
  labs(x = "Day of the year", y = "Predicted SST (°C)")

ggsave(paste0(home, "/figures/supp/BT_test_plots2.pdf" ), width = 17, height = 17, units = "cm")

# Plot with data
ggplot(filter(d, area == "BT" & year >= min(ndbt$year) & year <= max(ndbt$year)),
       aes(yday, temp, color = source)) + 
  geom_point(size = 0.001, alpha = 0.3) + 
  geom_line(data = ndbt, aes(yday, pred, color = source, linetype = "main model"), alpha = 0.7) +
  geom_line(data = ndbt2, aes(yday, pred, color = source, linetype = "BT specific model"), alpha = 0.7) + 
  scale_color_brewer(palette = "Dark2") + 
  facet_wrap(~year) + 
  labs(title = "Area = BT", color = "", linetype = "") + 
  guides(color = guide_legend(title.position = "top", title.hjust = 0.5)) + 
  theme_sleek(base_size = 8) +
  theme(legend.position = "bottom", 
        legend.direction = "horizontal",
        legend.spacing.y = unit(-0.3, "cm")) + 
  labs(x = "Day of the year", y = "Predicted SST (°C)")

ggsave(paste0(home, "/figures/supp/BT_test_plots3.pdf" ), width = 17, height = 17, units = "cm")

# Plot only new predictions, color by year facet by source?
p1 <- ggplot() + 
  geom_line(data = ndbt2, aes(yday, pred, color = year, group = year), alpha = 0.7, linewidth = 0.2) + 
  scale_color_viridis() + 
  facet_wrap(~source) + 
  labs(title = "Area = BT", color = "", linetype = "") + 
  guides(color = guide_legend(title.position = "top", title.hjust = 0.5)) + 
  theme_sleek(base_size = 8) +
  theme(legend.position = "bottom", 
        legend.direction = "horizontal",
        legend.spacing.y = unit(-0.3, "cm")) + 
  coord_cartesian(ylim = c(0, 25)) +
  labs(x = "Day of the year", y = "Predicted SST (°C)", title = "BT specific model")

p2 <- ggplot() + 
  geom_line(data = ndbt, aes(yday, pred, color = year, group = year), alpha = 0.7, linewidth = 0.2) + 
  scale_color_viridis() + 
  facet_wrap(~source) + 
  labs(title = "Area = BT", color = "", linetype = "") + 
  guides(color = guide_legend(title.position = "top", title.hjust = 0.5)) + 
  theme_sleek(base_size = 8) +
  theme(legend.position = "bottom", 
        legend.direction = "horizontal",
        legend.spacing.y = unit(-0.3, "cm")) + 
  coord_cartesian(ylim = c(0, 25)) +
  labs(x = "Day of the year", y = "Predicted SST (°C)", title = "main model")

p1 / p2

ggsave(paste0(home, "/figures/supp/BT_test_plots4.pdf" ), width = 17, height = 17, units = "cm")

# Now plot the annual means. Focus on 1980-2004, so that we don't need to worry about
# the Forsmark predictions in years before heating

ndbt_sum <- ndbt |> 
  filter(source == "logger") |> 
  group_by(year) |> 
  summarise(mean_pred = mean(pred)) |> 
  mutate(model = "main model")

ndbt2_sum <- ndbt2 |> 
  filter(source == "logger") |> 
  group_by(year) |> 
  summarise(mean_pred = mean(pred)) |> 
  mutate(model = "BT-specific model")

bt_sums <- bind_rows(ndbt_sum, ndbt2_sum)

ggplot(bt_sums, aes(year, mean_pred, color = model)) + 
  geom_line() +
  labs(x = "Year", y = "Model-predicted annual average temperature") + 
  scale_color_brewer(palette = "Accent")

ggsave(paste0(home, "/figures/supp/BT_test_plots5.pdf" ), width = 17, height = 17, units = "cm")

If we want to get uncertainty, we can use nsim instead; this simulates from the linear predictor using the inverse precision matrix, which is a fast way to get a distribution of samples from which we can take e.g. quantiles and means. However, it’s still slow, so the code below isn’t executed yet.

nd_sim <- data.frame(expand.grid(yday = seq(min(d$yday), max(d$yday), by = 1),
                                 area = unique(d$area),
                                 year = unique(d$year))) |>
  mutate(source = "logger") |>
  mutate(id = paste(year, area, sep = "_"),
         source_f = as.factor(source),
         year_f = as.factor(year))

# Trim!
nd_sim <- left_join(nd_sim, gdat, by = "area")

nd_sim <- nd_sim |>
  mutate(growth_dat = ifelse(year > min, "Y", "N")) |>
  filter(growth_dat == "Y") |>
  filter(yday %in% c(gs_min:gs_min)) |>
  mutate(area = as.factor(area))

# Predict!
nsim <- 500
m_pred_sims <- predict(m, newdata = nd_sim, nsim = nsim)

# Join sims with prediction data
nd_sim_long <- cbind(nd_sim, as.data.frame(m_pred_sims)) |>
    pivot_longer(c((ncol(nd_sim) + 1):(nsim + ncol(nd_sim))))

# Summarize sims over growing season
sum_pred_gs <- nd_sim_long |>
    ungroup() |>
    group_by(year, area) |>
    summarise(lwr = quantile(value, prob = 0.1),
              est = quantile(value, prob = 0.5),
              upr = quantile(value, prob = 0.9)) |>
    ungroup()

# In order to have a lower temperature in before-nuclear times (without any data to inform that), we can use the nearby areas..
sum_pred_gs <- preds |>
  mutate(keep = "Y",
         keep = ifelse(area == "BT" & year < 1980, "N", keep),
         keep = ifelse(area == "SI_HA" & year < 1972, "N", keep)) |>
  filter(keep == "Y")

sum_pred_gs_sub <- preds |>
  mutate(keep = "N",
         keep = ifelse(area == "FM" & year < 1980, "Y", keep), # use FM instead of BT
         keep = ifelse(area == "SI_EK" & year < 1972, "Y", keep)) |> # use SI_EK instead of SI_HA
  filter(keep == "Y")

# Now change the labels to BT and SI_EK...
sum_pred_gs_sub <- sum_pred_gs_sub |>
  mutate(area = ifelse(area == "FM", "BT", "SI_HA"))

# Bind rows and plot only the temperature series we will use for growth modelling
sum_pred_gs <- bind_rows(sum_pred_gs, sum_pred_gs_sub) |> select(-keep, -type)

order <- sum_pred_gs |>
  group_by(area) |>
  summarise(mean_temp = mean(temp)) |>
  arrange(desc(mean_temp))